Tutorial Part IV - Movement - Solution Exercise 2

Exercise 4.2

Three animals have been collared and tracked, and we want to know their home range size and whether the animals actually avoid agriculture and prefer forest.
There is telemetry data available from three animals (DieTiere_32633.shp) in the folder data-raw. Additionally, land cover (Landnutzung_3035.shp) and rivers (linearegew_32633.shp) are available in the geo-raw folder. Land use has a column ‘NS1’ with coded info on 1 (urban areas), 2 (agriculture), 3(forest), 4(swamps), 5 (water bodies).

After loading libraries, the data and making a plot, please do the following:
a) Calculate their home ranges.
b) Buffer the home range with the mean displacement per step of all three animals
c) Clip the home ranges with the underlying landscape and calculate the percentage
of land use classes 2/ agriculture and 3/ forest within the home range (hint: use st_intersection)
d) In which land use are the animals actually located? Extract land use underneath each telemetry location (hint: st_join). Do the animals actually avoid agriculture?

Get started

# Open libraries
## general basics
library(here)
library(ggplot2)

## the basic packages for spatial analyses and viz
library(sf)
library(terra)
library(tmap)

## for telemetry and home range analyses
library(move)
library(adehabitatHR) 
library(sp)

Define workspace and folder structure

work_wd <- here()
data_raw_wd <- here('data-raw')
geo_raw_wd  <- here('data-raw/geo-raw')
output_wd   <- here('output')

Load data

# if you are not sure about the file names, check the folder content with
# list.files(data_raw_wd)
animals <- st_read(dsn=data_raw_wd,layer="DieTiere_32633")  
## Reading layer `DieTiere_32633' from data source 
##   `C:\Users\kramer\PopDynIZW Dropbox\Steph Kramer\_GitHub\Course4_MoveQ\data-raw' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 388 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 390397.5 ymin: 5862042 xmax: 396748.2 ymax: 5868366
## Projected CRS: WGS 84 / UTM zone 33N
landcov <- st_read(layer="Landnutzung_3035", dsn=geo_raw_wd)    
## Reading layer `Landnutzung_3035' from data source 
##   `C:\Users\kramer\PopDynIZW Dropbox\Steph Kramer\_GitHub\Course4_MoveQ\data-raw\geo-raw' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 182 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 4539362 ymin: 3304861 xmax: 4566087 ymax: 3324860
## Projected CRS: ETRS89-extended / LAEA Europe
rivers  <- st_read(layer="linearegew_32633", dsn=geo_raw_wd)   
## Reading layer `linearegew_32633' from data source 
##   `C:\Users\kramer\PopDynIZW Dropbox\Steph Kramer\_GitHub\Course4_MoveQ\data-raw\geo-raw' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2834 features and 13 fields
## Geometry type: MULTILINESTRING
## Dimension:     XYZ
## Bounding box:  xmin: 381989.2 ymin: 5852103 xmax: 408856.6 ymax: 5871804
## z_range:       zmin: 0 zmax: 0
## Projected CRS: WGS 84 / UTM zone 33N

Define or transform CRS

lc_32633 <- st_transform(landcov, 32633)

Make a plot of the data

#Create a color code for the 5 land use classes
mycol <- c('red','antiquewhite','chartreuse4','cornsilk4','blue')

## TODO - argh, Punkte und Linien erscheinen nicht im html....

plot(lc_32633[,'NS1'], col=mycol[lc_32633$NS1],border='transparent')
plot(animals[,1],col=animals$ID,add=TRUE) 
plot(rivers[,1], col= 'blue', lwd=0.1,add=T)

#Create a color code for the 5 land use classes
mycol <- c('red','antiquewhite','chartreuse4','cornsilk4','blue')

# TODO - plot mit colorcode mycol funktioniert nicht....

ggplot(lc_32633) + 
  geom_sf(data=lc_32633[,'NS1'], aes(fill  = mycol[NS1]), color='transparent') +
  geom_sf(data=animals) + 
  geom_sf(data=rivers, color='blue') +
  theme_classic()

#Create a color code for the 5 land use classes
mycol <- c('red','antiquewhite','chartreuse4','cornsilk4','blue')
myid  <- c('yellow','pink','red')

tmap_mode(mode = "view")
tm_shape(shp = lc_32633[,'NS1']) + 
  tm_polygons(col = "NS1", palette = mycol) +
tm_shape(shp = animals[,'ID']) + 
  tm_dots(col = "ID", palette = myid) +
tm_shape(shp = rivers[,1]) + 
  tm_lines(col = 'blue')

a) & b)

First , we retrieve the steplength. For this, transform the sf-Object into a move-Object.

anim_coords <- as.data.frame(st_coordinates(animals)) #returns a list
animals$X <- anim_coords$X
animals$Y <- anim_coords$Y
animals_ng <- st_drop_geometry(animals)

animals_move <- move(x=animals_ng$X, y=animals_ng$Y,
                     time = as.POSIXct(animals_ng$MYDATUM, 
                                       format = "%Y-%m-%d 8:0:0", tz="CET"),
                     proj=CRS("+init=epsg:32633"),
                     data=animals_ng, animal = animals_ng$ID, sensor='GPS')
plot(animals_move, type='l') 

# since we have 3 animals, a list with 3 elements is returned
steplength_per_animal <- distance(animals_move)
mean_stepl_p_anim <- lapply(steplength_per_animal,mean)
(mean(unlist(mean_stepl_p_anim)))
## [1] 303.0478

The mean step length is approx. 300 m. Now use this value to buffer the MCPs.

# for using mcp() in {adehabitatHR}, convert to SpatialPointsDataframe of {sp}
animals_sp <-as(animals, "Spatial")
mcp <- mcp(animals_sp[,'ID'], percent=95) 

# the buffer command needs however an sf object
mcp_sf  <- as(mcp, "sf")
mcp_buf <- st_buffer(mcp_sf, 300)

# the plot is correct
plot(mcp)                         
plot(animals_sp, add=TRUE, col= animals_sp$ID) 
plot(mcp_buf[,1], add=TRUE)
plot(animals_sp, add=TRUE, col= animals_sp$ID) 
plot(mcp, add=TRUE)

# save buffer
st_write(obj = mcp_buf, 
         dsn = output_wd, 
         layer = 'mcp_buffer_32633',
         driver = 'ESRI Shapefile', 
         delete_layer = TRUE)
## Deleting layer `mcp_buffer_32633' using driver `ESRI Shapefile'
## Writing layer `mcp_buffer_32633' to data source 
##   `C:/Users/kramer/PopDynIZW Dropbox/Steph Kramer/_GitHub/Course4_MoveQ/output' using driver `ESRI Shapefile'
## Writing 3 features with 2 fields and geometry type Polygon.

c)

# check the total area of the buffer in km2 per animal:
buf_area <- st_area(mcp_buf)/1e6

# clip with land use
mcp_buf_lc <- st_intersection(lc_32633, mcp_buf) 
# have a look at the object:
mcp_buf_lc
## Simple feature collection with 45 features and 14 fields
## Geometry type: GEOMETRY
## Dimension:     XY
## Bounding box:  xmin: 390158.2 ymin: 5861742 xmax: 397048.2 ymax: 5868666
## Projected CRS: WGS 84 / UTM zone 33N
## First 10 features:
##        AREA PERIMETER CLCGES_BES CLCGES_B_1  NS BRD SZENE   SPHEROID   TKNR NS1
## 2  33397200  79888.10      23163      23162 312   1 19323 Krassovsky N33111   3
## 5  11562800  44001.70      23265      23264 211   1 19323 Krassovsky N33111   2
## 15 18208500  73976.60      23764      23763 231   1 19323 Krassovsky N33111   2
## 16  8568800  37160.80      23788      23787 231   1 19323 Krassovsky N33111   2
## 28  5142930  13701.10      24028      24027 312   1 19323 Krassovsky N33111   3
## 35  1227710   4956.00      24314      24313 243   1 19323 Krassovsky N33111   2
## 36  1448920  11624.70      24348      24347 311   1 19323 Krassovsky N33111   3
## 41   654073   3560.28      24431      24430 211   1 19323 Krassovsky N33111   2
## 42  3534710  11425.00      24454      24453 311   1 19323 Krassovsky N33111   3
## 44 46474200 125933.00      24480      24479 211   1 19323 Krassovsky N33111   2
##    NS2 NS3 id     area                       geometry
## 2   31 312  1 1444.323 MULTIPOLYGON (((394196.9 58...
## 5   21 211  1 1444.323 POLYGON ((392577.5 5866547,...
## 15  23 231  1 1444.323 POLYGON ((391413.6 5866716,...
## 16  23 231  1 1444.323 POLYGON ((395139.3 5866196,...
## 28  31 312  1 1444.323 POLYGON ((390833 5866283, 3...
## 35  24 243  1 1444.323 POLYGON ((392292.4 5866469,...
## 36  31 311  1 1444.323 POLYGON ((395809.1 5865790,...
## 41  21 211  1 1444.323 POLYGON ((394729.2 5865833,...
## 42  31 311  1 1444.323 POLYGON ((390395.2 5865319,...
## 44  21 211  1 1444.323 MULTIPOLYGON (((392247.2 58...
# take care, you have to recalculate the area!
mcp_buf_lc$areakm2 <- st_area(mcp_buf_lc)/1e6

# now, for each animal id, summarize the area per NS1-category 2 and 3:
agri   <- tapply(mcp_buf_lc$areakm2[mcp_buf_lc$NS1 == 2],
               mcp_buf_lc$id[mcp_buf_lc$NS1 == 2],sum)
forest <- tapply(mcp_buf_lc$areakm2[mcp_buf_lc$NS1 == 3],
                 mcp_buf_lc$id[mcp_buf_lc$NS1 == 3],sum)

# percentages
round((agri/buf_area*100),0)
## Units: [1/m^2]
##  1  2  3 
## 71 48 85
round((forest/buf_area*100),0)
## Units: [1/m^2]
##  1  2  3 
## 29 52 14

d)

spatjoin <- st_join(animals,lc_32633)
myres <- table(spatjoin$NS1,spatjoin$ID)
myres
##    
##       1   2   3
##   2  48  25  95
##   3 147  64   9
# now divide that by the total number of locations per animal
myn <- table(animals$ID)

round((myres[1,]/myn*100),0)  # agriculture
## 
##  1  2  3 
## 25 28 91
round((myres[2,]/myn*100),0)  # forest
## 
##  1  2  3 
## 75 72  9
# compare the values with the percentages above

For example, animal 1 had 71% of agriculture and 29% of forest in its HR, however, only 25% of the locations were actually in agriculture, but 75% of its locations were in forest. We could conclude that animal 1 avoids agriculture.

END

Session Info 
Sys.time()
## [1] "2022-03-07 02:33:24 CET"
sessionInfo()
## R version 4.1.2 (2021-11-01)
## Platform: i386-w64-mingw32/i386 (32-bit)
## Running under: Windows 10 x64 (build 19043)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_Germany.1252  LC_CTYPE=English_Germany.1252   
## [3] LC_MONETARY=English_Germany.1252 LC_NUMERIC=C                    
## [5] LC_TIME=C                       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] adehabitatHR_0.4.19 adehabitatLT_0.3.25 CircStats_0.2-6    
##  [4] boot_1.3-28         MASS_7.3-54         adehabitatMA_0.3.14
##  [7] ade4_1.7-18         deldir_1.0-6        move_4.1.6         
## [10] rgdal_1.5-28        raster_3.5-11       sp_1.4-6           
## [13] geosphere_1.5-14    tmap_3.3-2          terra_1.4-22       
## [16] sf_1.0-5            ggplot2_3.3.5       here_1.0.1         
## 
## loaded via a namespace (and not attached):
##  [1] RColorBrewer_1.1-2      httr_1.4.2              rprojroot_2.0.2        
##  [4] tools_4.1.2             bslib_0.3.1             utf8_1.2.2             
##  [7] R6_2.5.1                KernSmooth_2.23-20      DBI_1.1.1              
## [10] colorspace_2.0-2        withr_2.4.3             tidyselect_1.1.1       
## [13] leaflet_2.0.4.1         compiler_4.1.2          leafem_0.1.6           
## [16] textshaping_0.3.6       xml2_1.3.3              bookdown_0.24          
## [19] sass_0.4.0              scales_1.1.1            classInt_0.4-3         
## [22] proxy_0.4-26            systemfonts_1.0.3       stringr_1.4.0          
## [25] digest_0.6.29           rmarkdown_2.11          base64enc_0.1-3        
## [28] dichromat_2.0-0         pkgconfig_2.0.3         htmltools_0.5.2        
## [31] highr_0.9               fastmap_1.1.0           htmlwidgets_1.5.4      
## [34] rlang_0.4.12            farver_2.1.0            jquerylib_0.1.4        
## [37] generics_0.1.1          jsonlite_1.7.3          crosstalk_1.2.0        
## [40] dplyr_1.0.7             magrittr_2.0.1          s2_1.0.7               
## [43] Rcpp_1.0.8              munsell_0.5.0           fansi_1.0.2            
## [46] abind_1.4-5             lifecycle_1.0.1         stringi_1.7.6          
## [49] leafsync_0.1.0          yaml_2.2.1              tmaptools_3.1-1        
## [52] grid_4.1.2              parallel_4.1.2          crayon_1.4.2           
## [55] lattice_0.20-45         stars_0.5-5             knitr_1.36             
## [58] pillar_1.6.4            codetools_0.2-18        wk_0.5.0               
## [61] XML_3.99-0.8            glue_1.6.0              evaluate_0.14          
## [64] leaflet.providers_1.9.0 png_0.1-7               vctrs_0.3.8            
## [67] rmdformats_1.0.3        gtable_0.3.0            purrr_0.3.4            
## [70] assertthat_0.2.1        cachem_1.0.6            xfun_0.28              
## [73] lwgeom_0.2-8            e1071_1.7-9             ragg_1.2.1             
## [76] class_7.3-19            viridisLite_0.4.0       tibble_3.1.6           
## [79] memoise_2.0.1           units_0.7-2             ellipsis_0.3.2